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ABSTRACT 

We demonstrate the capability of Principal Component Analysis (PCA) as applied by Heyer 
& Schloerb (1997) to extract the statistics of turbulent interstellar velocity fields as measured by 
the energy spectrum, E{k) ~ . Turbulent velocity and density fields with known statistics 
are generated from fBm simulations. These fields are translated to the observational domain, 
T(x,y,v), considering the excitation of molecular rotational energy levels and radiative transfer. 
Using PCA and the characterization of velocity and spatial scales from the eigenvectors and 
eigenimages respectively, a relationship is identified which describes the magnitude of line profile 
differences, 5v and the scale, L, over which these differences occur, Sv ~ L". From a series of 
models with varying values of /3, we find, a — 0.33/3 — 0.05 for 1 < < 3. This provides the 
basic calibration between the intrinsic velocity field statistics to observational measures and a 
diagnostic for turbulent flows in the interstellar medium. We also investigate the effects of noise, 
line opacity, and finite resolution on these results. 

Subject headings: hydrodynamics — turbulence — ISM: kinematics and dynamics — ISM: clouds line: 
profiles — methods: statistical 

Introduction 



The evolution of the cold, dense interstellar medium and the formation of stars are regulated by the 
kinetic energy content of turbulent gas flows. Such chaotic motions play a critical role in the equilibrium of 
the gas as these provide a countering pressure to self gravity and the pressure of the external medium. In the 
absence of external pressure fronts, the formation of stars occurs exclusively in dense, localized regions where 
the turbulent energy is sufficiently reduced to enable gravitational collapse. Therefore, an understanding of 
the dense interstellar medium and the formation of stars requires an accurate description of turbulent gas 
flow including energy injection into and dissipation from the gas system and the distribution of turbulent 
kinetic energy over varying spatial scales. 

Turbulent flows within the molecular interstellar medium are recognized as observed supersonic line widths 
and morphological complexity (Scalo 1990; Falgarone Phillips, & Walker 1991). While these observations 
denote the presence of turbulence in the interstellar medium, there is little quantitative information as to the 
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nature of the turbulent flows. Ideally, one would like to recover the S-dimcnsional velocity field, Vz{x,y, z), 
for the line of sight component, z, from a set of spectroscopic observations. In practice, the reconstruction 
of the true velocity field from observations is not possible given the chaotic, non-deterministic nature of 
turbulent flows and the non-linear transformation of the field onto the spectroscopic axis of the observer. 
Therefore, it is essential to construct statistical descriptions of turbulence (Dickman 1985). The primary 
statistic used to distinguish turbulent models is the energy spectrum, E(k), which describes the degree of 
coherence of the velocity field over spatial scales. In D dimensions, the energy spectrum is the angular 
integral of the power spectrum, 

E{k) = j dflPik) oc< P{\k\) >Q k^-^ (1.1) 

We refer to the spectral index j3 in connection with E(k), 

E{k) ~ Ifcl-'^ (1.2) 
The power spectrum then has the form (in D dimensions) : 

P{k) ~ (1.3) 

For a dissipationless cascade of energy through an incompressible fluid, /? = 5/3 and the power spectrum, 

P(k), decreases in 3 dimensions as k"^-'^/'' (Kolmogorov 1941). A 2 or 1 dimension cut through the 3 
dimensional structure will also have an energy spectrum with the same index j3 but with a different slope of 
the power spectrum. 

The mean square turbulent velocity, vi, at size I, is determined from the spectrum of velocity fluctuations 
over scales smaller than 

/>00 /'OO 

~ / E{k)dk (X / k-'^dk cx l/^'^ (1.4) 

vf describes the variance of the 3 dimensional velocity field, v(x,y,z) within a volume l^. The root mean 
square turbulent velocity, < vf >^/^, is then, 

< vf >V2oc = p (1.5) 

where 7 = (/3 — l)/2. It is not equivalent to an observed width of a spectral line which depends upon the 

velocity, density, temperature, and abundance which vary along the line of sight. However, given a velocity 
field described by an energy spectrum with /? > 0, one might expect some correlation of measured velocities, 
5v, with projected spatial scale, L, of the form, 

6v ~ L" (1.6) 

where a refers to the measured spectral index in contrast to the intrinsic index 7 in equation (1.5). Conven- 
tionally, 6v has been derived from measured line width or the displacement of centroid velocities of spectral 
lines. 

Coarse limits to values of P appropriate for the molecular interstellar medium can be inferred from the 
early models of observed CO line profiles. To account for the broad velocity widths and the ability of the 
high opacity CO lines to detect the warm molecular gas located deep within the Orion cloud, Goldreich & 
Kwan (1974) modeled line profiles with a purely systematic velocity field replicating collapsing motions. For 
such smooth, systematic velocity fields, /? » 3. However, it is unlikely that all molecular clouds undergo 
global gravitational collapse since this predicts too large a star formation rate in the Galaxy (Zuckerman & 
Evans 1974). Other systematic velocity fields such as rotation are not observed over the extended scales of 
molecular clouds (Arquilla & Goldsmith 1984). These results suggest /3 must be smaller than 3. A lower 
limit to p is provided by the inability of purely random velocity fields {(3 = -2 in 3D) to reproduce observed 
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line profiles (Leung & Lizst 1976; Dickman 1985). Given these considerations, an approximate range to the 
spectral index of the energy spectrum is 3/2 < (3 < 5/2. 

More direct measures of P rely on the identification of a size-line width relationship and the assumption 
that variations along the spectroscopic axis accurately reflect variations within the intrinsic velocity field. 
Using a set of inhomogeneous data from the literature, Larson (1981) determined a relationship between the 
global velocity dispersion, tT„, and size, L, for a sample of 46 molecular regions in the solar neighborhood, 
ay ~ L^-^^. Using a more homogeneous sample of observations, Solomon etal (1987) identified a relationship 
with an index of 0.5. However, the applicability of these observed relationships to the energy spectrum 
associated with interstellar turbulence (equation (1.5)) is weak unless one considers the targeted clouds as 
part of singular, global fluid within the Galaxy (see Dickman 1985). Given the wide range of distances 
and location of the clouds within different spiral arms, such an assumption is clearly invalid. Rather, these 
respective relationships derived from a sample of unrelated clouds provide a statement of self gravitational 
equilibrium if one assumes that the mean column density of clouds is approximately constant. That is, there 
is a sufficient amount of mass distributed within the cloud to bind the internal turbulent motions. 

Several investigations have decomposed spectroscopic data cubes of ^^CO or C^*0 emission from 
targeted molecular clouds into discrete units or "clumps" based on some preassigned definition (Carr 1987; 
Stutzki & Gusten 1990; Williams, de Geus, & Blitz 1994). For each clump, a global velocity dispersion 
and size are determined. However, the resultant size-line width relationships show a large scatter and little 
significant correlation. This may reflect that the localization of emission within a spectroscopic data cube 
may not uniquely originate from a discrete clump localized in space. Also, the process of identifying clumps 
limits the resolution to the size of the smallest object. 

Scalo (1984) and Kleiner & Dickman (1985) proposed the use of structure and autocorrelation functions 
of centroid velocity images derived from spectroscopic data cubes. This assumes that the 2 dimensional 
projection retains the statistical characteristics of the 3 dimensional velocity field. These correlation methods 
have been applied by Miesch & Bally (1994) to a sample of 12 molecular regions imaged in ^^CO emission. 
The structure functions are described by a power law with a spectral index of 0.86 which corresponds to 
a=0.43 as in equation (1.6). However, these correlation methods are limited. For many regions in the 
molecular ISM, observed line profiles exhibit non-gaussian shapes for which the centroid velocity is not an 
accurate approximation. Also, the centroid velocity is not well defined near the edges of the clouds where 
the signal is weak. Finally, the relationship between the measured value of a and the intrinsic velocity field 
statistics, 7 or /3, has not yet been established. 

Recently, Heyer & Schloerb (1997, hereafter HS97) motivated the use the multivariate technique of Princi- 
pal Component Analysis (PGA) to decompose spectral line imaging observations of molecular clouds. From 
the analysis, they identify a relationship which describes the magnitude of velocity differences of line profiles 
and the spatial scale at which these differences occur within the image. For several molecular clouds, they 
find a ~ 0.4 — 0.5. However, it has not yet been demonstrated that this analysis, or the other methods 
described above, are sensitive to varying field statistics (/?) nor has the relationship between a and /3 been 
established. Given the complex transformation of the velocity field to the observed spectral axis, one can 
not assume that a = 7 = (/? — l)/2. 

In this paper, we quantitatively evaluate the ability of PGA to characterize intrinsic velocity fields from 
spectral line observations. Guided by theoretical and observational results, 3 dimensional models of velocity 
and density fields with known statistical properties are generated. Observations of these fields are simulated, 
aecoimting for non-LTE excitation and radiative transfer. The simulated observations are analyzed with PGA 
under a variety of physical conditions and observational limitations. We demonstrate that the method is 
indeed sensitive to velocity fields with varying energy spectra and derive the appropriate calibration between 
the observed quantity a and /3. In a companion paper (Brunt & Heyer 2000), we use these results to derive 
the energy spectra for molecular regions in the outer Galaxy. 
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2. Construction of Model Turbulent Clouds 



2.1. Velocity Field Generation 

To evaluate the sensitivity of the method outUned by HS97 to varying velocity field statistics, we construct 
3 dimensional model velocity fields with prescribed energy spectra. Fields which exhibit varying degrees of 
spatial correlation can be generated from "colored noise" or fractional Brownian motion (iBm) simulations 
(Voss 1988; Stutzki et al. 1998). These simulations are not solutions to the fiuid equations but rather, 
are expedient numerical tools with well defined spatial statistics. The resultant fields are oversimplifications 
of real interstellar velocity fields in which one expects shock structures, intermittency and other complex 
features. 

We use the simplest version of an fBm vector velocity field, in which the components are generated 
independently. For this type of field there is a roughly equal amount of power between the divergence free 
(solenoidal) and the curl free (compressible) components. Our model velocity fields are thus somewhat more 
extreme (in terms of the energy fraction in the compressible modes) than what may be expected for 'typical' 
turbulent fields which tend to have their energy spectra dominated by the divergence- free component. The 
effects, if any, of differing energy fractions in each mode will be the subject of a future paper. 

The model velocity field generation is a straightforward process which exploits the Fourier space imple- 
mentation of fBm in three dimensions. The three steps to generate the v(x,y,z) field are: (1) distribute 
delta-correlated (Gaussian) noise on a cubic grid of N x N x N pixels. In this study, we use N=128; (2) 
Fourier transform this field to the frequency domain, filter the Fourier space field to the desired power law 
energy spectrum of spectral index P, and transform back to the real space domain; (3) Normalize the field 
to the desired variance, a^. In this study, all velocity fields are normalized such that = 1. Only one 
component of the vector velocity field, the observable line-of-sight component 'tu^{x.y, z), is generated. For 
our analysis, we consider purely isotropic fields. This means that the longitudinal gradients Svz/Sz are 
statistically equivalent to the transverse gradients {Svz/Sx and Svz/Sy). 

Ten realizations of v(x,y,z) fields are generated for each of several values of /3 (1.0, 1.5, 2.0 ,2.5, 3.0, 4.0). 
The input noise for each realization has a unique random number seed. Representative one-dimensional cuts 
through fields with varying energy spectra generated by this method are displayed in Figure 1. It is clear 
that as /3 increases, the velocity field becomes smoother and similar to a linear gradient. 

2.2. Real-Space Velocity Field Statistics for Model Fields 

To ensure that the generated fields have the desired real-space statistics (i.e. 7), we made direct (real- 
space) measurements on the generated fields. Of course, it would be possible to construct full-resolution 
autocorrelation functions or structure functions but without the use of the Fourier Transform, these would 
be lengthy calculations. We would like to avoid the use of Fourier space in checking the fields, since these 
were created within the Fourier domain. There is also the requirement that the correct normalization for the 
measurements be found. To check the real-space statistics, we calculate a variant of the standard (second 
order) structure function S{t), 

S{t) =< {v{r) - v{r + T)f > (2.1) 
where the angle brackets indicate averaging over all points separated by lag r. We seek the desired dependence 

5(t) oc r^T (2.2) 

S{t) can be usefully approximated via coarse-grained versions of the original field (Brunt 1999 and references 
therein). Specifically, the field v is partitioned into an ensemble of cubical boxes of L pixels on a side, over 
which V is averaged. This procedure is carried out for dyadic values of L (1, 2, 4, 8, 16, 32, 64), where L=l 
corresponds to the original field. For each coarsened version of the field, the mean square difference along 
the cardinal directions between adjacent coarse pixels is calculated, providing an estimate of S(t = L). The 
coarse structure function 

Sc{L)=<ivL{r)-VL{r + L)f > (2.4) 
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where vi, is a coarse-grained field. Prom these measurements the exponent 7 is obtained via the proportion- 
ahty 

Sc{L) oc (2.5) 

The ensemble coarse structure function, < Sc{L) >, is formed by equal weight averaging of all Sc(L) 
measurements at each L for each set of ten fields at each (3. These ensemble measurements are shown in 
Figure 2, The deficit at L=64 is partly caused by the periodicity of the fields. A power law is fit to < Sc{L) > 
for L < 32 to determine an ensemble value, "Ie, for each /?. Also, a mean value, < 7 > is derived from the 
average of all 7 values obtained individually at each (3. These results are summarized in Table 1. It is clear 
that the calculated values of 7 are equivalent to the expected value (/3 — l)/2, derived from the input value 
of (3. In Figure 3, we show the variation of 7^ with the intrinsic value of (3. For velocity fields with /3 > 3, we 
find that 7 asymptotically approaches unity. This is due to the rapid decrease of the power spectrum such 
that only the |k| = 1 components are significant. Therefore, the field is similar to a smooth linear gradient. 
In a recent study, Myers & Gammie (1999) derive a similar conclusion. The lower values of 7^: and < 7 > 
relative to the expected values for the larger (3 values are likely caused by the periodicity requirement, which 
suppresses the out-of-band (large-scale) gradients, and also partly due to the gradients which result from 
stochastic fluctuations of the power in the first two wavenumbers. That is, an excess of power at the first 
wavenumber cannot make the field more smooth (higher 7), but a deficit of power at this wavenumber can 
make the field less smooth (lower 7). A better approximation to "absolutely smooth" fields requires the 
extension to /3 = 4. 

2.3. Density Field Generation 

We consider two forms of the density field. The basic calibration is established with a uniform density field, 
no(a;, y, 2)=constant. However, to provide a more realistic observation, an inhomogeneous density field is 
also considered. Since the spatial variation of density c;orrc;sponds to varying excitation of the trac;c; molecule, 
these inhomogeneous models are used to to investigate the effects of non-uniform sampling of the velocity 
field. Recent hydrodynamical models for an isothermal gas show a lognormal density probability density 
function, PDF (Padoan et al 1997: Passot & Vazquez- Semadeni 1998). Additionally, the power spectrum of 
the density distribution, in the simulations of Padoan et al (1997), is consistent with a power-law form, of 
spectral slope (in the power spectrum) of ~ 2.6 and similar to values determined from observations (Stutzki 
et al. 1998). 

To generate a field that has a lognormal PDF and a power-law power spectrum, we implement a modified 
version of the fBm method described in Section 2.1. An £Bm field with /3 = 1 is exponentiated which 
produces a resultant field with the desired statistics (Schertzer & Lovejoy 1987; Pecknold et al. 1993). The 
(3=1 fBm field serves as Inn -|- constant, where the normalizing constant, given the dispersion in the field, 
determines the mean density. The normalization of the field to a mean density of rio = 1 is obtained by 
setting < Inn >= — (Tinn, where dinn is the global standard deviation of Inn. Note that dinn is independent 
of no. Like the iBm velocity fields, the density field is periodic across the cube boundaries. Therefore, 
the density field is inspected and, if necessary, re-centered to ensure that no major structures lie across the 
boundary. We emphasize that there is no spatial correlation between the density and velocity fields as one 
might expect for real clouds in which density fluctuations result from the localized divergence of the velocity 
fleld. 

2.4. Physical Properties of the Models 

In order to calculate realistic opacities, we assign a physical scale to the simulation. For all realizations 
the linear size of the cube is 20pc. A pixel or cell is 20pc/128=0.156 pc. The mean density for the uniform 
density field is 100 cm~^. For the log normal filds, the mean density varies from 10^ to 10^ cm~^. The gas 
is assumed to be isothermal with a kinetic temperature, T^ equal to 20 K. An approximate measure of the 
global gravitational equilibrium of the models is given by the ratio of kinetic (~ 3cr^ Por^) to gravitational 
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(~ GpqT^) energies 

'AalJGpar^ - 1.2(100/ < uh^ >) (2.6) 

where we have taken the radius, r= 10 pc. The models do not significantly deviate from gravitational 
equilibrium although the equilibrium state of the model clouds is not critical to these results. 

2.5. Observations of the Models 

Given the density and velocity fields and the kinetic temperature, the brightness temperature is determined 
from the excitation of the chosen gas probe and the radiative transfer of the emission toward the observer's 
line of sight. The excitation of the line, parameterized by the excitation temperature, Tex, and line center 
opacity. To, at the position (x,y,z), are determined from a non LTE calculation which accounts for local 
radiative trapping (Scoville & Solomon 1974). This depends on the local kinetic temperature, density, and 
Nfl-jX/Av, where N//^ is the molecular hydrogen column density, X is the abundance ratio to molecular 
hydrogen, and Av is the local line width. For these simulations, we generate synthetic line profiles of ^'^CO 
J=l-0 emission, including excitation calculations up to J=5. Constant values of ^^CO abundance relative to 
H2 (1.25 X 10~^) and a kinetic temperature (20 K) are assumed. A pixel scale of (20 pc/128) = 0.156 pc 
gives the H2 column density within a cell at position (x,y,z) as 

NH2{x,y,z) = ZmxlQ^^ n{x,y,z) {2Q/I2d,) cm'^ (2.7) 

The local line width is determined from the quadrature sum of the thermal velocity for the molecule, Cth, 
and the gradient of the velocity field at that position 

I^v{x, y, z) = ^I^~^^?^EmA^2 (2.8) 

Figure 4 shows the variations of excitation temperature and line center opacity with density and Nj^^X/Av 
for the assumed kinetic temperature of 20 K. Once the excitation temperature and line center opacity are 
determined for each point in the cloud, the emergent intensity at velocity, v, is calculated for each position 

(x,y) 

T{x,y,v) = -f- V(i5,(Te,) -B,(2.7))(l -e-^°*(^(^'^'^)'^''))e-^(^'^'^'") (2.9) 

where B,^(T) is the Planck function, v is the line frequency, (^(v(x,y,z),Av) is the Gaussian profile function 
normalized to unity and 

z 

r{x, y, z,v)=Y^ r{x, y, C)Hv{x, y, C, M) (2-10) 

C=o 

is the total foreground opacity to a physical depth z. The model cloud is placed at a distance of 1 kpc and 
is observed with a top-hat telescope beam response which covers 1 pixel or an effective angular resolution of 
32.4". The spectral resolution of the observation is 0.05 kms~^. Examples of the synthetic observations can 
be found in Brunt (1999). 

3. Analysis 

The formal goal of PGA is to determine the set of orthogonal axes such that the data, when projected 
onto these axes, maximizes the variance. In practice, for spectroscopic imaging observations, it effectively 
identifies differences in line profiles as these contribute to the variance of the data cube. Such line profile 
differences can arise from dynamical processes within the gas and can be distinguished from the instrumental 
noise of the data. Therefore, PGA provides a powerful method to identify kinematic variability within the 
target object. In this section, we describe the PGA technique and subsequent structural analysis to extract 
the observationally accessible exponent a from spectroscopic imaging observations. 
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3.1. Principal Component Analysis Method 

A spectroscopic imaging observation is comprised of an ensemble of n = nxXriy spectra each with p 
spectroscopic channels. We write the resultant data cube as T{xi, yi, vj) = T{ri, vj) = Tij where rj = (a;,, t/j) 
denotes the spatial coordinate of the ith spectrum. The covariance matrix Sjk is 

1 " 

Sjk = - '^TijTik (3.1) 

The set of eigenvectors, uij and eigenvalues are determined from the solution of the eigenvalue equation for 
the covariance matrix, 

Su = Xu (3.2) 

The eigenvalue, A;, corresponds to the amount of variance projected onto its corresponding eigenvector, 
Ui. Therefore, the eigenvalues (and corresponding eigenvectors) are reordered from largest to smallest. In 
this implementation of PC A, we do not explicitly subtract the mean value, < Tj >= Y^Tij/n from each 
channel image. These are effectively removed within the first, {I = 1), principal component. Note that 
no spatial information is contained within the eigenvectors, since the ordering of the spatial part of Tij is 
arbitrary. The eigenvectors Uij are purely spectroscopic and trace the (ordered) sources of variance in the 
ensemble of observed line profiles. Also, the eigenvectors are orthogonal, which ensures that the information 
contained within each component is independent information. Similar decompositions can be obtained by 
other orthogonal function sets such as wavelets or spherical harmonics. The advantage of PCA over these 
other basis sets is that the orthogonal vectors are determined by the data and not predetermined. 

To spatially isolate the sources of variance contained with the l*^ component, an eigenimage, /'(r,), is 
constructed from the projected values of the data, Tij, onto the eigenvector, uij, 

p 

l\ri) = J2TijUij (3.3) 
j=i 

It is equivalent to a velocity-integration of T(x,y,v) that is weighted at each channel by the value of the 
appropriate eigenvector. Since PCA is a linear decomposition, the random noise of the input data can 
readily be propagated. The variance of projected values due to noise is 

amjf = f;^ulamjf (3.4) 

.7=1 

Assuming, that cr{Tij) is constant over the bandpass of the spectrum and recalling that u is orthonormal 
{J^'^'ij = 1)' the expression reduces to 

aiin) = a{Tij) (3.5) 

Thus, the variance of values projected onto the component is equal to variance of the input data. More 
importantly, one can readily distinguish values within the eigenimage from those due to instrumental noise. 

3.2. Characteristic Scale Measurements 

For a given component, I, the eigenvectors and eigenimages describe the velocities and positions at which 
the measured line profiles are different with respect to the noise level. Characteristic velocity differences 
and the spatial scales over which these differences occur are derived from e-folding lengths of the normal- 
ized autocorrelation functions (ACFs) of the eigenvectors and eigenimages respectively (Heyer & Schloerb 
1997). The raw autocorrelation functions, Cy{dv), C|(r), are calculated directly from the eigenvectors and 
eigenimages respectively, 

C{r{dv) =< {u\v)u\v + dv) > (3.6a) 



drin) 
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C| (r) =< (/' (r)/' (r + r) > (3.66) 

However, in this form, the 2 dimensional autocorrelation fimction contains an additive component, Cn{t), 
due to instrumental noise. In Section 3.3, we demonstrate that the noise subtracted autocorrelation function, 
Cjq(t) = Cj{t) — Cn{t) where Cn{t) is the autocorrelation function of noise values of the data. 

The characteristic velocity scale, 6vi, is determined from the velocity lag at which Cy((5?;j)/Cy(0) = e~^. 
For the 2 dimensional ACF of the eigenimage, one must additionally correct for the effects of finite resolution 
observations upon the zero lag (see Section 3.4). We determine the spatial correlation lengths along the 
cardinal directions, dxi,5yi, from the noise corrected ACF, 

= e-^ (3.86) 
and assign the biased characteristic spatial scale. Lib to the quadrature sum, 



Lib = yj5x1 + 5y'^ (3.9) 



As discussed in Section 3.4, the true characteristic scale, Li, is derived from the Lis with subpixel corrections 
to account for finite resolution of the observations. Since the cigcnimages tend to be isotropic, the ACF 
profiles along the cardinal directions are a reasonable approximation to the variation of the ACF along any 
arbitrary angle. The characteristic velocity and spatial scales are statistically well defined with respect to 
the noise as those represent mean values over the full spectrum and field respectively. A measurement error, 
ugy, to the value &vi is given by one half the velocity resolution of the spectrometer. For the spatial scale, 
aL is determined from the quadrature sum of the spatial resolution and the degree of anisotropy. \Sxi — 5yi\, 
in the spatial ACF. The PCA velocity statistic, a, is then obtained from the set of Jw, L pairs which are 
larger than the respective spectral and pixel resolution limits and form the power law relationship 

where the unsubscripted quantities refer to the ensemble {5vi,Li). Note that in the retrieval of each Sv — L 
pair, PCA incorporates information from the entire data cube such that each pair is well defined. 



3.3. Instrumental Noise Effects 

All real observations contain a noise contribution. For any analysis, it is critical to evaluate how the noise 
may aff'ect the result. For the ACF of the eigenvector, the characteristic velocity scale measurements are 
largely imaffected by instrumental noise due to the limited number of channels (p<128). However, for large, 
2 dimensional images, the contributions to the ACF due to noise can be significant and must be removed. 
The observed eigenimage is the sum of a noiseless eigenimage, Io{r), and a noise contribution, N{r), 

f{r)=ll,{r) + Nir) (3.10) 

The unnormalized spatial ACF, (C|(t)) is 

C|(r) =< (4(r) + A^(r))(4(r + r) + N{r + r)) > (3.11) 

C\{t) = CUt) + C%{t)+ < 4(r)iV(r + r) + ^(r + r)N{r) > (3.12) 

where Cjg is the ACF of the noiseless eigenimage and C^v is the ACF of the noise. The term in brackets 
measures the correlation between signal and noise. We assume that the signal and noise are uncorrelated 
and treat this term as zero. Therefore, 

C\o{r) = C\{r)-&^{r) (3.13) 



8 



For Gaussian noise, 

rl (-r\ - f '^N for T = 

<~^N[r) - I Q for T 7^ 

We exploit the noise propagation properties of PCA and obtain cr^ from several high I eigenimages, which 

are free of signal. A more general procedure can be followed to remove the noise contributions to the 
autocorrelation function. The noise ACF, Cjyf(r), can be calculated directly from either channel images 
within the data cube or high I eigenimages which contain no signal. The noise contribution is removed at 
each value of r. This method works well for non-Gaussian noise fields such as those obtained by on the fly 
mapping or reference sharing observations (Brunt & Heyer 2000). 

The noise subtracted ACF, Cjg is well determined even in the presence of relatively high noise levels. If 
this noise subtraction is not carried out, then Cj contains a noise spike at the origin, such that the spatial 
correlation lengths are underestimated. Since the noise contribution relative to the signal is larger for higher 
I eigenimages, this means that the scale lengths are underestimated in such a way as to lessen the slope 
a of the of the log((5?;)-log(iy) relationship. HS97 did not make this correction for the noise contribution. 
Therefore, the derived values of a obtained by HS97 should be considered as lower limits. 

3.4. Correction for Finite Resolution Effects 

It is important that the analysis results are independent of the resolution at which the observations are 
obtained. In this Section we discuss the systematic effects within the ACFs due to finite resolution. We 
demonstrate in Section 4.4 that the spectroscopic ACFs and the derived velocity scales are not compromised 
by changing the spectroscopic resolution at which the observations are obtained. However, the eigenimage 
ACFs are resolution-dependent and require a sub-pixel correction to ensure that the spatial scale measure- 
ments are not systematically compromised. In brief, for finite resolution observations, the unnormalized 
ACF at zero lag, C|q(G), is imderestimated. However, by assuming, or estimating, a shape to the ACF in 
order to extrapolate to zero lag, one can gauge the amount that the scale length is overestimated and provide 
a correction to the derived biased scale lengths, Lib- The alternative option is to fit a functional form with 
a "scale" parameter to the ACF. We have found that making small corrections to scales obtained from a 
biased ACF provides the more stable measurement. 

A detailed derivation of the corrections for various functional forms of the ACF is given by Brunt (1999). 
In brief, the measured value of C/o(0) which is used to normalize the ACF is 

Cio{0) « C{eLpi,) (3.14) 

where C(r) is the true ACF, Lpi^ is the pixel size, and e < 1, is a number that depends on the shape of 
the ACF and the shape of the telescope beam. For eigenimages derived from observations with a top hat 
beam profile, e ~ 0.5. This underestimation with respect to an ideal ACF applies to the full resolution 
(128^) eigenimages, since a simulation or observation at any finite resolution necessarily involves omission of 
variance (i.e. contributions to C/o) that would be present in a higher resolution measurement. Fortunately, 
the degree of underestimation and the associated scale corrections for the 128^ eigenimages are small. 
The overestimated scale measurement Lb obtained from an underestimated ACF, do, satisfies 

^ = e.p(-l) (3.15) 

For a functional form of the true ACF, C 

C{t) = exp{-{T/LY) (3.16) 

where L is the true scale length of the ACF and k parameterizes the shape of the ACF, the biased scale 
length is is related to the true scale length L via 

L^{L%- (eLpi,)«)V« (3.17) 
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with e « 0.5. For n = 1 (an exponential ACF), the correction to all measured scales is —eLpi^ ~ —Lpix/2. 
For a Gaussian beam —eLpix « —0.75Lpix (Brunt 1999). The form of the ACFs for each observation is 
estimated from a few low I eigenimages. The corrections are not sensitively dependent on small errors in 
estimating k. The application of these corrections to the measured scale lengths ensures consistency of the 
results for varying resolutions. 

4. Results 

Initially, the basic calibration between the observed exponent a and the intrinsic field statistics, param- 
eterized by /3, is established for uniform density fields. Subsequently, we evaluate the 5v — L relationship 
with respect to the uniform density results using more realistic distributions of interstellar material and 
observational conditions including instrumental noise and varying resolution. 

4.1. Uniform Density Results 

Ten realizations of the velocity field are created for each value of /3 ranging from 1.0 to 4.0 according to 
the prescription given in Section 2.1. The real space statistics of these fields are chocked and reported in 
Section 2.2. Synthetic line profiles of ^^CO J=l-0 emission are calculated at each position, (x,y), assuming 
a uniform density (n//^ = 100 cm~^) and uniform temperature (Tfc = 20K) via the LVG radiative transfer 
method described in Section 2.5. The optical depths in these observations are small and the excitation 
is subthermal and uniform throughout the field which corresponds to an almost perfectly-sampled velocity 
field. The resultant data cubes are decomposed using PCA and velocity and spatial scales are estimated from 
the eigenvectors and eigenimages respectively. These scales define a relationship between the magnitude of 
velocity differences and the spatial scales at which these differences occur. The relationship is fit to a power 
law parameterized by the exponent a. 

Figure 5 shows 5v — L relationships for 6 velocity fields with varying values of (3. The dotted lines mark the 
spectroscopic and spatial pixel sizes. Uncorrected spatial scales Lb and spectroscopic scales 6v are restricted 
to > 2 pixels and > 1 spectroscopic channel respectively (see Section 3.4). The lower (3 simulations are 
characterized by a fewer number of retrieved 5v — L pairs above the resolution limits since most of the 
variance is generated at small scales. 

In Table 1 and Figure 6 we show the relationship between the mean value of a determined from the 
ensemble of realizations and the input /3. The plotted error bars reflect the standard deviation of values 
from this mean. These resiilts show that there is a monotonic relationship between the observed and intrinsic 
exponents for /3 < 3 and the method, as described in Section 3, is sensitive to velocity fields with varying 
statistics. A least squares fit to the a, (3 values for 1 < /3 < 2.5 gives 

a= (0.33 ± 0.04)/3 - 0.05 ±0.08 (4.1) 

This was obtained with a as the dependent variable. For values of /3 > 3, the derived values of a show 
a similar transition to smoothness that was seen in the intrinsic 7 measurements described in Section 2.2. 
Equation 4.1 provides the basic calibration between an observable measure of kinematic variability within a 
data cube and the intrinsic velocity field statistics. The relationship between a and (3 remains an empirical 
measurement. It surely depends upon the dimensionality of the fields since a recent calibration of this 
technique to 2 dimensional fields reveals a modified relationship, 020 ~ fhu/^- Giv(ni the extreme complexity 
of transforming a velocity field onto a spectroscopic axis and the identification of velocity gradients by line 
profile differences, an analytical description of equation 4.1 is not readily determined. 

4.2. Lognormal Density Field Results 

The results of the previous Section are an important first step in gaining insight into how PCA retrieves 
the intrinsic velocity field information. However, a uniform density field is clearly unphysical. To evaluate 
the method with a more realistic density distribution, lognormal density fields (see Section 2.3) with mean 
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values of 10^, 10"^, and 10^ cm~'^ and velocity fields with /3 = 1, 2, 3 are observed. The increase of the mean 
density corresponds to an increase in gas column density and line center optical depth within a given cell. 
The primary effect of inhomogeneous density fields is to mask the velocity field from the observer as the low 
density regions do not sufficiently excite the molecular gas tracer into emission. Figure 7 shows the Sv — L 
relationships obtained by PCA from these observations for these varying conditions. The results obtained 
from the uniform density observations with the identical velocity fields are included for comparison. The 
resultant a measurements are tabulated in Table 2. 

For the /3 > 2, the derived values of a are similar to the uniform density exponents although there is a 
constant offset in the logarithm. In addition, the number of significant components at the largest spatial 
scales increases with increasing density and optical depth. This is surely due to radiative trapping within 
a given cell which maintains an elevated excitation temperature in low density regions in which the trace 
molecule would otherwise be insufficiently excited. This suggests that a molecular gas tracer with moderate 
to high optical depths, such as ^^CO, provides a more complete sampling of the velocity field. The /3 = 1 
observations tend to overestimate a relative to the uniform density exponent. In this case, the presence of 
high frequency velocity fluctuations places many cells along a given line of sight at a common velocity (see 
Figure 1). Therefore, for optical depths larger than unity, the back side of the cloud is hidden from the 
observer and the field is undersampled. From these results, we conclude that the velocity field statistics with 
inhomogeneous density fields can be reliably recovered for /3 > 2 under a wide range of optical depths but 
are slightly overestimated for more shallow energy spectra /3 ~ 1. 

4.3. Effects of Instrumental Noise 

Any real observation of the ISM necessarily contains instrumental noise. In this Section, we evaluate the 
effects of instrumental noise upon the analysis by adding a noise field, N(r,v), to the simulated observations 
of the lognormal density field (< >= 10^ cm~"^) with values of /? = 1,2,3. The input noise is Gaussian 
with variance a%. We define C as a gauge of the signal to noise of a data cube. 



where ct^ = (Tq + '^n ^^le variance of the data cube with noise and ctq is the variance of the same in a 
noiseless observation. Since PCA generally consolidates all spectroscopic channels with signal within the 
first principal component and noise in the higher components, ( can be calculated from the I = 1 and I = p 
eigenvalues, 



In Figure 8, we show the Sv — L relationships derived from simulations with varying signal to noise 
(C = 4.0,2.0,1.0). Also shown are the relationships derived with ( = oo (dotted line). The resultant a 
measurements are shown in Table 3. The primary effect of noise is to reduce the number of significant 
components that are detected above the resolution limits as the variance generated by the velocity field 
which is contained in the higher components are overwhelmed by the variance due to noise. This tends to 
induce systematic effects at large /? rather than an increased scatter in the retrieved Sv — L pairs. For most 
spectral line imaging observations, ^ > 4 for which the results are largely unaffected by noise (Brunt & Heyer 
2000). 

4.4. Resolution Effects 

In Section 3.6, we described the necessary corrections to the retrieved spatial scales which minimize 
the effects of finite resolution observations. In this section, we demonstrate the need and utility for such 
corrections using the simulated observations with varying resolution. The original simulated observations 




(4.2) 




(4.3) 
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are comprised of 128 x 128 spatial pixels. We degrade the resolution of the simulated observation (/3 = 2 
and lognormal density field with < n(i?2) >= 10'* cm~^) by spatial block-averaging. This is a reduction 
in the number of pixels and is equivalent to observing the simulation at a further distance with a top hat 
beam. We define A/ as the number of resolution elements across the image. Block averaged versions of 
this simulated observation are calculated for A/ = 64,32, 16. Further reduction in resolution results in only 
one detected measurement above the resolution limits. The measurements obtained from these observations, 
for both uncorrected scales and scales corrected for an exponential ACF (k = 1) are shown in Figure 9. 
These plots demonstrate the saturation of retrieved pairs at the limit L = (1 — exp{—l))Lpix, expected 
for a linear interpolation between unity and zero in the normalized ACF. This feature is characteristic of 
uncorrelated noise with a true scale length of zero. The uncorrected measurements clearly show steepening 
due to a uniform scale overestimation as the resolution limit, Lpi^, is approached. It is this steepening 
that compromises the estimation of a. When the corrections are applied to the retrieved spatial scales, the 
derived values of a are consistent between the varying resolutions. Note that if no corrections are made, then 
the measured value of a represents an upper limit, obtained with the assumption that there is no sub-pixel 
inhomogeneity. Measurements on real data will always contain this ambiguity, which can only be resolved 
by increasing the resolution. 

It was stated in Section 3.4 that the velocity scales 5v are unaffected by changing the spectroscopic 
resolution. To demonstrate this statement, we coarse-grained the velocity axis of the same observation. In 
analogy to the above discussion of spatial scales, we define A„ as the ratio of the total range in velocity 
spanned by the spectrometer and the velocity resolution. For the original simulated observation, A„ = 128. 
The spatial resolution is fixed at A/ = 128. The results of the PCA retrievals for these fields are shown 
in Figure 10. Also shown are scale measurements for a coarse-grained field with A/ = 64, and A„ = 32. 
The spatial scales have been corrected for finite spatial resolution. No corrections have been applied to the 
retrieved velocity scales. These results shows that the velocity scale measurements with degraded spectral 
resolution are equivalent to those derived from the full resolution observations. 

5. Summary 

We have investigated the ability of PCA to recover intrinsic turbulent velocity field statistics from spectral 
line imaging observations of the molecular ISM, and established a preliminary calibration of the method. 

The calibration, subject to the limitations of the modeling, identifies a monotonic relationship between 
the intrinsic f3 values (or equivalently, 7) and observed exponents (a) that characterize the velocity field 
statistics. 

In addition to the simplistic nature of the input density and velocity fields, including the assumption 
of statistical isotropy and the periodicity feature, the major limitations are the lack of mutual consistency 
between density and velocity, and the assumption of local excitation in the radiative transfer calculations. 
However, the simulated lognormal observations are in good visual correspondence with real molecular ISM 
observations. The retrievals from the low density observations are in close agreement with those from the 
high density observations. This suggests that the accuracy of the method depends upon the degree to which 
one samples the velocity field, rather than the details of the radiative transfer. For molecular clouds in 
particular, highly saturated observations may be more reliable (due to more homogeneous sampling and a 
greater number of retrieved measurements) and this suggests that ^^CO is likely to provide the most accurate 
values of a. This method could also be applied to HI 21cm emission which is excited over a wide range of 
physical conditions. 

The sampling effects due to noise, which uniformly masks sources of variance, do not affect the fits for 
a, until very high noise levels are reached. As calibrated, the PCA method is resolution-independent but 
limited in general by necessary assumptions of sub-pixel homogeneity of the emission field. Independent of 
our modeled calibration, the robustness of measured a values provide quantitative constraints to simulated 
observations of velocity and density fields from hydrodynamical calculations. We do stress, however, that 
the calibration of a to /3 achieved in this paper is strictly only valid for the fBm model velocity fields. More 
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realistic turbulent fields may not adhere as closely to the calibration, and in future we intend to include 
supercomputer MHD simulations in the PCA calibration. 

Empirical measurements of a from real data of course do not depend on the calibration to /3 (or 7), 
and such measurements will be reported by Brunt & Heyer (2000) for an ensemble of fields taken from the 
FCRAO Outer Galaxy Survey. Using the preliminary calibration established here, we obtain constraints on 
the intrinsic velocity field statistics of the molecular ISM. 

This work was supported by NSF grant AST 97-25951 to the Five College Radio Astronomy Observatory. 
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Fig. 1. — Examples of one- dimensional cuts through the 3 dimensional model velocity fields for 1 < /3 < 4. 
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Fig. 2. — Ensemble coarse structure function measurements of the model velocity fields. Major tick marks 
correspond to powers of ten. 



16 




13 

Fig. 3. — Measured relationship between (3 and 7 derived from the coarse structure function. 
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Fig. 4. — Dependence of T,,.,-, (left) and Tq (right) on the local (pixel) density, column density and velocity 
gradient. The To wedge is logioarithmic. 
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Fig. 5. — Representative results of the PCA retrieval for the uniform density observations. The dotted lines 
mark the spatial pixel size and spectroscopic channel width. 
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Fig. 6. — Basic calibration of the observable a. exponent to the intrinsic exponent /3. 
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Fig. 7. PCA mcasiircincnts from the observation of the lognormal density field. The dashed lines denote 
the fits obtained from the observation of the same velocity fields with a uniform density distribution. The 
solid lines are the fits to 5v-L pairs. 
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Fig. 8. — PCA measurements from observations of the lognormal density field with varying degrees of 
gaussian distributed noise. 
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Fig. 9. — Spatial resolution effects on the PCA measurements. The As = 128 measurements are shown as 
filled circles, and the labeled lower resolution versions as open circles. The dotted lines mark the "full" 
resolution limits and the solid line the lower resolution limits. 
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Fig. 10. — Spectroscopic resolution effects on the PCA measurements. Ttie = 128 measurements are 
shown as filled circles, and the labeled lower resolution versions as open circles. The dotted lines mark the 
"full" resolution limits and the solid line the lower resolution limits. 
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Table 1: Uniform Density Field 





^= ^ 

' 2 


IE 


< 7 > 


a 


1.0 


0.0 


0.00±0.03 


-0.02±0.04 


0.30±0.06 


1.5 


0.25 


0.26±0.02 


0.23±0.04 


0.45±0.03 


2.0 


0.50 


0.50±0.01 


0.48±G.05 


0.62±0.02 


2.5 


0.75 


0.72±0.01 


0.73±0.04 


0.79±0.04 


3.0 


1.00 


0.85±0.02 


0.86±0.05 


0.92±0.04 


4.0 


1.50 


0.93±0.02 


0.94±0.01 


1.02±0.06 



Table 2: Lognormal Density Fields 



/3 




a 






n = lO^cm ^ 


n = 10'*cm~^ 


n = lO^cm"^ 


1.0 


0.40±0.01 


0.47±0.01 


0.44±0.01 


2.0 


0.57±0.03 


0.61±0.02 


0.64±0.01 


3.0 


0.92±0.03 


0.97±0.03 


0.92±0.03 



Table 3: Lognormal Density Field with Noise 
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a 










C = 00 


C = 4.0 


C = 2.0 


C = i.o 


1.0 


0.44±0.01 


0.45±0.01 


0.41±0.03 


0.48±0.01 


2.0 


0.64±0.02 


0.64±0.03 


0.64±0.02 


0.57±0.02 


3.0 


0.92±0.03 


0.88±0.02 


0.81±0.01 


0.63±0.04 
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